Potential energy landscape and long time dynamics in a simple model glass 



L. Angelani^, G. Parisi^, G. Ruocco^ and G. Viliani"^, 
^ Universitd di L'Aquila and Istituto Nazionale di Fisica della Materia, 1-67100, L'Aquila, Italy. 
^ Universitd di Roma La Sapienza and Istituto Nazionale di Fisica Nucleare, 1-00185, Roma, Italy. 
Universitd di Trento and Istituto Nazionale di Fisica della Materia, 1-38050, Povo, Trento, Italy. 

(February 1, 2008) 

We analyze the properties of a Lennard- Jones system at the level of the potential energy landscape. 
After an exhaustive investigation of the topological features of the landscape of the systems, ob- 
tained studying small size sample, we describe the dynamics of the systems in the multi-dimensional 
configurational space by a simple model. This consider the configurational space as a connected 
network of minima where the dynamics proceeds by jumps described by an appropriate master 
equation. Using this model we are able to reproduce the long time dynamics and the low tem- 
perature regime. We investigate both the equilibrium regime and the off-equilibrium one, finding 
those typical glassy behavior usually observed in the experiments such as: i) stretched exponential 
relaxation, ii) temperature-dependent stretching parameter, lii) breakdown of the Stokes-Einstein 
relation, and iv) appearance of a critical temperature below which one observes deviation from the 
fluctuation-dissipation relation as consequence of the lack of equilibrium in the system. 



PACS Numbers : 61.20.Lc, 64.70.Pf, 82.20.Wt 
I. INTRODUCTION 

The landscape paradigm is a very useful point of 
view for the study of glassy systems. The detailed anal- 
ysis of the free energy or potential energy surfaces allows 
us to get insight into the rich phenomenology exhibited 
by glass-forming liquids in the supercooled phase, around 
the glass transition region and in the low temperature 
glassy regime. While the investigation of the free energy 
surface is a very hard task starting from the microscopic 
description of the system, since the landscape details are 
very strongly temperature dependent, the description of 
the potential energy landscape is a much more tractable 
problem, and is a good starting point to investigate prop- 
erties at not-too-high temperatures. The trajectory of 
the representative point of the system in the configura- 
tional phase space can be viewed as a path in the mul- 
tidimensional potential energy surface. The dynamics is 
strongly influenced by the topography of the landscape: 
local minima, barriers height, attraction basins and fur- 
ther topological features. In recent years it has been show 
that the details of the potential energy surface are of 
great importance in determining the properties of many 
systems exhibiting glassy behavior, like glass-forming liq- 
uids, protein folding, atomic cluster or evolutionary bio- 
logical models. 

In this paper we numerically investigate the low tem- 
perature dynamical properties of a simple system, i.e. a 
monoatomic Lennard- Jones system, through an analy- 
sis of its multidimensional potential energy surface and 
a simple model for the loe temperature dynamics. In 
the first part (Section II), by studying small sample size, 
we give a comprehensive description of the potential en- 
ergy landscape of the systems: minima, barriers, reaction 
paths, saddle points, and determine statistical distribu- 
tions and cross-correlations among the analyzed quanti- 



ties. In the second part (Section III) we use this infor- 
mation to set up a simple model for the study of the 
long-time relaxation dynamics of the system; the model 
consists of a connected network of potential energy min- 
ima with a jump dynamics described by an appropriate 
master equation. This allows us to get information on the 
behavior of the system for times long enough that a di- 
rect Molecular Dynamics (MD) simulation is not feasible 
Q . After a static (thermodynamical) test of the model, 
we determine the dynamical equilibrium properties and 
also the off-equilibrium ones, and discuss the results. In 
the last Section we report the conclusions. 

II. POTENTIAL ENERGY LANDSCAPE 

We numerically investigate the topology of the poten- 
tial energy hypersurface of a Lennard- Jones 6—12 system 
of N interacting particles in a cubic box with periodic 
boundary conditions. The pair potential is: 



with r the Euclidean distance between two particles. The 
physical parameters a and e are choose to describe an Ar- 
gon system: a = 0.3405 nm and e/Ks = 125.2 K {Kb 
is the Boltzmann constant). In order to have an exhaus- 
tive description of the energy landscape, we investigate 
small size systems, with a number of particles N < 30. 
As we shall see later, such small systems can nevertheless 
exhibit complex enough behavior. 

Due to the small size of the system, the range of in- 
teractions among the particles is of the same order as 
the box length. It is then much more appropriate to use 
a multi-image method instead of the usual minimum im- 
age method H , in which each particle interacts only with 
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the nearest image (generated by the periodic boundary 
conditions) of aU other particles. In our case many im- 
ages contribute, and it is necessary to consider them in 
the calculation of the interactions. The choice we made 
introduces an angular dependence in the pair potential 
which, however, is negligible with respect to the radial 
one. The simulated density is p = 4.2 • 10~^ mol/cm^, 
which was obtained by Demichelis et al. Q as the small- 
est value at which a Lennard- Jones system of = 864 
Argon particles is stable in the glassy state after a rapid 
quench from high temperature. 

A. Minima 

As a first step we search for the potential energy min- 
ima. We used a modified conjugate gradient method 
starting from high temperature configurations obtained 
during a MD simulation. In this way we find the so 
called "inherent configurations", corresponding to local 
minima. New minima are identified by their potential 
energy values. We analyzed systems with A^ = 11 29 
and, for each N, we stopped the search when the rate at 
which new low energy minima are found (for example in 
lowest first third of the full energy range) is smaller than 
a given number (about 10~^). In this way we were able 
to obtain a good classification of the low energy minima. 
The number of found minima, Af, does not show a clear 
and well defined dependence on A, contrary to the case 
of small clusters Q , so that it is not possible to give an 
estimate of the coefficient a in the expected exponential 
growth: Af oc exp(Q;A^) (for clusters a ~ 1). We ob- 
serve indeed for some values of A^ a strong tendency of 
the system to fall always in the same minima. In this 
case an exhaustive research of the inherent structures is 
a very hard numerical task, this explains the not clear 
dependence found. 

Once most minima are classified, we analyze their fea- 
tures: energy, static structure factor, curvature and dis- 
tances. 

1. Energy 

The first important information about an inherent 
configuration is its potential energy value. From the 
energy distribution it is sometimes possible to recog- 
nize the crystalline-like configurations. Usually one finds 
an evident gap between the lowest energy minima, the 
crystalline- like ones, and the other ones with higher en- 
ergies, as can be seen in Fig. 1 where we report the 
distribution of minima energies per particle for A^ = 29. 
In other cases the situation is not so clear and it is use- 
ful to use a different method to characterize the nature of 
an inherent configuration; the most useful quantity in de- 
termining the spatial order structure of a minimum is its 



static structure factor, which allows (imperfect) crystal- 
like configuration also of high energy to be identified. 

2. Static structure factor 

In order to classify the spatial distribution of particle 
in a given minimum of the potential energy we use the 
static structure factor: 

Siq)^j^\J2''M^'l-rj)\' ■ (2) 

j 

Due to the finite size of the system the allowed q vectors 
are of the form q — (2-k/L) n, with n — {n,j;,ny,nz) an 
integer vector. We define the quantity 

5(g) = - ^(9^' (3) 

qe{q,Ag} 

where J2q£{q Aq} is ^ o^^'" '^i vectors with mod- 
ulus within q ± Ag. For a pure crystalline configuration 
of A^ particles S{q) consists of Bragg peaks, and its value 
at the peaks is Smax = N. For amorphous configurations 
S{q) does not present a well defined peak structure and 
the highest value is Smax ~ 2 -f- 3, obtained for a q value 
of the order of the inverse mean distance of two near par- 
ticles. For small sized systems, like those analyzed here, 
there are intermediate situations and we use the criterion 
Smax < A/2 to determine the amorphous nature of an 
inherent structure. 

3. Curvature 

Another important property of a minimum is its overall 
curvature c, defined as the determinant of the Hessian $ 
of the potential energy function <I>: 

c = det($"), (4) 

The eigenvalues of the Hessian matrix are proportional 
to the squared vibrational eigenfrequencies. In the inset 
of Fig. 1 we show the distribution of 

^ = Logio (c/m)i/2 , (5) 

where m is the mass of the particles (we use the Argon's 
value m = 40 amu). The quantity 7 is thus proportional 
to the sum of the logarithms of the frequencies of normal 
vibrational modes 

a 

with a = 1, Af — 3 (the three zero frequencies corre- 
sponding to rigid translations have been eliminated from 
the sum). As can be seen in Fig. 1, the highest 7 values 
correspond to the minima with lowest energy, i.e. the 
crystalline-like minima which are narrower and deeper 
than the other packing structures (see also Sec. II-C). 
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4- Stress tensor 

For each minimum we determined the ofF-diagonal part 
of the microscopic stress tensor: 



i<j 



•^ij Uij 



(7) 



where Xij , yij , Zij are the components of rij ~ ri — fj ; 
these quantities wiU be useful in determining the shear 
viscosity. The form of the stress tensor we use is the 
q ~ extrapolation of the g-dependent expression |^] 



•exp(ig • n) , 



with 



1 — exp(ig • r) 
iq ■ r 



(8) 



(9) 



The index i and j refer to particles, while a and (3 label 
the spatial axes x, y, z. In (|^) we have omitted the kinetic 
term, not well defined in an inherent configuration; our 
hypothesis is that also this truncated form is able to well 
describe the relaxation processes. 



5. Distances 

We now turn to the relationships among the minima; 
in particular we determined the mutual distances in the 
3iV-dimensional configuration space; in this regard it is 
important to take into account all the symmetry opera- 



tions of the problem. Indicating with = (f^ , . 



' N 



the 3N coordinates of the particles in the minimum a, 
we define the distances dab between minima a e b: 



lab = min(|r„ -rf,|) 

I .H.,7T 



(10) 



where the minimization is made with respect to the con- 
tinuous translations (T), discrete rotations and reflec- 
tions (i?), and permutations (tt) of the particles. 

The minimization over the continuous translations (T) 
is done by putting, sequentially, each particle of a in the 
same place as each one of b. The minimization over the 
rotations and reflections (R) is carried out by consider- 
ing the 48 symmetry operations of the cubic group; the 
minimization over particle permutation (tt) is apparently 



an hard task to solve, as one should consider all the N\ 
possible configurations in a direct calculation, but this is 
not the case actually. The problem is of polynomial type, 
i.e. the time needed to find the optimum solution grows 
as a polynomial function of the size, N, of the system 
(on the other hand non-polynomial problems require an 
exponential computational time, i.e. the travel salesman 
problem). The optimization problem we need to solve 
is a bipartite matching problem, which can be done in 
very short computational time by using an appropriate 
algorithm. 



B. Barriers 

A very important topological quantity in determin- 
ing the dynamical behavior of the system is the energy 
barrier experienced by the system in traveling from one 
minimum to another. At first sight it might seem that 
it is accurate enough to evaluate the barrier along the 
straight path joining two minima in the BA^-dimensional 
space. However, as evidenced by Demichelis et al. [|], 
this produces in most cases much higher barriers than 
an evaluation along the least action path, indicating that 
the straight path approximation is often not good. We 
have then determined the least action path for each pair 
of minima a and b, which is defined as the path that 
minimizes the action functional 



S[£] = / dsy/2m[<i>{r{s)) 
Ji 



$0 



(11) 



where € is a generic path between the minima, s is the 
curvilinear coordinate and $o = niin{$a, ^b}- This func- 
tional problem is simplified by dividing the path into a 
finite number of intervals (typically we use n = 16 inter- 
vals) and by minimizing the action function with respect 
to the extreme of the n segments constrained to move 
in hyperplanes perpendicular to the straight path. The 
highest potential-energy value along the least action path 
determines the barrier height and identifies the saddle 
point. 

We have only analyzed the system with < 17, due 
to the very long computational times needed in the cases 
with A^ > 17 where there are too many pairs of minima 
to take into account. We report the results for the largest 
system analyzed, A^ = 17 with TV = 38 minima. In Fig. 2 
we show an example of the potential energy profile along 
the straight path (dashed line) and along the least action 
path (full line) between two minima. It is evident that 
the energy barrier is significantly lower in the latter case, 
although the two paths are not very distant in configu- 
ration space. Similar results are obtained for the other 
analyzed paths. Sometimes it happens that two minima 
are not directly connected, in the sense that the least ac- 
tion path joining them crosses a third minimum, and a 
non trivial connectivity among the minima emerges (in 
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the analyzed system we find that each minimum on an 
average is directly connected to 20 other minima). For 
each saddle point along the least action path we deter- 
mine the main properties: energy, curvature, down "fre- 
quency" . In Fig. 3 is reported the energy distribution 
of the barriers A$fcar. • The curvature is defined as the 
absolute value of the determinant of the Hessian of the 
potential energy evaluated at the saddle: 



R 



I det $ 



sad. I 



Csad. = I det $ 



sad. I 7 



(12) 



In the inset of Fig. 3 we show the distribution of the 
quantity 



(13) 



The down "frequency" uJsad. is defined as the square root 
of the absolute value of the down curvature along the 
least action path: 



(14) 



where v is the tangent vector to the least action path at 
the saddle point. 



C. Correlations 

In order to have a full statistical description of the 
potential energy landscape, it is useful to investigate, be- 
sides the distributions of the different quantities, also 
their cross-correlations. We have then determined the 
linear correlation coefficient 



-x){y^ -y) 



(15) 



for all the measured quantities x and y. In Table I are 
reported the values obtained, together with the log — log 
correlations, in order to evidence possible power laws. 
In the first column we report the number of particles of 
the system analyzed (only for TV = 17 we have deter- 
mined the least action paths and all the related quan- 
tities). It emerges that, substantially, energy difference 
and distance among the minima are not correlated, in- 
dicating that the topological structure of the inherent 
configurations is not energy-correlated. A weak corre- 
lation is observed between energy and curvature at sta- 
tionary points (minima and saddle points). In Fig. 4 
we show the cross-correlation between the energies of the 
minima and their curvatures. An interesting correlation 
is observed between barrier energies A^har. and distances 
among minima dai, (Fig. 5), with a nearly linear correla- 
tion in double log scale (line in the figure). 

We conclude the analysis of the energy landscape by 
determining the entropy ratio R, defined as the ratio be- 
tween the curvature of saddle points and that of the re- 
lated minima 



det C 



(16) 



This quantity gives a quantitative measure of the abil- 
ity of the system in finding the right path to reach an- 
other minimum. If i? ^ 1 there are no entropic hin- 
drances, while if i? 1 these effects become relevant, 
as the narrowness (higher value of the curvature) at the 
saddle makes the least action path toward that specific 
minimum unfavorable with respect to other escape route. 
For all the minimum-saddle-minimum triplets we have 
evaluated R; the majority of the values is in the range 
10~^-^10, in qualitative agreement with the results found 
in Lennard- Jones clusters |q|. 



III. MODEL FOR THE DYNAMICS 

The investigation of the properties of glass-former liq- 
uids at the level of the energy landscape, allows us to 
introduce some approximations in the dynamics of the 
system. We define a simplified model which is able to 
capture the long time behavior of the system, and which 
consists of a connected network of potential energy min- 
ima with a jump dynamics among them described by an 
appropriate master equation 

The basic idea is quite simple. A glass is represented 
by a configurational point confined in a very small region 
of the accessible phase space and in the zero tempera- 
ture limit (neglecting quantum effects) all the atoms are 
frozen in well defined positions, corresponding to some 
mechanically metastable state. When the temperature 
is raised jumps among different mechanically stable posi- 
tions become possible. At finite and not-too-high temper- 
ature we assume the dynamically relevant processes are 
the following: a short time dynamics dominates by small 
vibrations around stable positions (this dynamics can be 
described within the harmonic approximation by diago- 
nalizing the dynamical matrix), and a long time dynam- 
ics consisting of collective (involving many atoms) jumps 
among different stable positions. The main hypothesis 
we make is that there is a substantially clear separation 
of time scales between the two dynamical processes. This 
characterization of the dynamics is a good approximation 
at not-too-high temperature. By increasing the tempera- 
ture, anharmonic effects become relevant to the vibration 
around the local minimum and, moreover, a clear time 
scale separation between fast vibrational and slow jumps 
dynamics is no longer possible. In a recent work Q the 
validity of this hypothesis has been verified in a Lcnnard- 
Jones binary mixture, giving a very good agreement with 
a direct MD investigation. To sum up, our model, which 
is expected to capture the physics of the system at low 
temperature, is based on two main hypotheses: 
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1. clear cut difference between vibrational dynamics 
at short time and dynamics of collective jumps at 
long time; 

2. description of the long time dynamics through a 
master equation, with the transition rates that de- 
pend on the topological properties of the potential 
surface. 

The main advantages of the introduced model with 
respect to the usual MD computations are: 

• we can avoid in a simple way the crystallization 
process, that always takes place in one component 
LJ systems, as we do not consider the crystalline 
minima in setting up the network. 

• we can study in a direct way the low temperature 
properties, where usually the very long relaxation 
times require very long computational time. In MD 
the computational times are proportional to the 
physical times, while in the model introduced here 
the computational times are those needed to find 
the eigenvalues and eigenvectors of the transition 
matrix, independent of temperature; 

• it is possible to evidence the relationships between 
the energy landscape and the behavior of the sys- 
tem. 

To be more specific, the model is a connected network 
of potential energy minima and the master equation gov- 
erning the jumps dynamics is: 



Wab pI = Wba 



(20) 



dpa 

dt 



{t;b,to)=Y, WacPc{t;b,to) 



(17) 



where Pa(t;b,tQ) is the probability that the system is at 
minimum a at time t, if it was at minimum b at time 
to. The off-diagonal elements of the matrix W are the 
transition rates. The diagonal elements are fixed by the 
condition 



(18) 



In order to obtain an asymptotic behavior that repro- 
duces the right Boltzmann weight the occupation proba- 
bility must satisfy: 

lim pa{t;b,to)=pl = ^(detC)-!/^ exp(-/3$,) , 

(19) 

{Z is such that J2a Pa — ^ ^^'^ the pre-exponential factor 
follows from the harmonic vibration in each minimum), 
and the transition matrix W must satisfy the detailed 
balance relation: 



The solution of the master equation is easily expressed 
in terms of eigenvalues A„ and eigenvectors ai""* {n = 
1, M, with M matrix dimension) of W: 

Pait;b,to) = ip',)-' "i"' "i"' exp[A„(i-io)] • 



(21) 

In the numerical calculus it is more convenient to express 
the solutions in terms of the eigenvectors of a new sym- 
metric matrix Wab = Wab [Pb/P^aY^"^ (whose eigenvalues 
coincide with those of W): 

Pa[t-bM) = [pVpIY'^ E 4"' exp[A„(t-io)] . 

n 

(22) 

The model is well defined once we give an appropriate 
form to the transition matrix W . In order to deter- 
mine the transition rates let us analyze the problem of 
escape from a metastable state; a useful point of view for 
systems with many degree of freedom is the description 
in term of a few relevant coordinates. This reduction 
is possible whenever there are few reaction coordinates 
with characteristic evolutionary times longer than those 
of the other degrees of freedom, which act as effective 
terms on the relevant coordinates, i.e. like noise and vis- 
cous terms. We suppose this is the case of our system 
whenever the temperature is not too high (the analysis 
of reaction paths made by Demichelis et al. H] supports 
this hypothesis). Handling the problem as a Markovian- 
Brownian d-dimensional motion in the overdamped fric- 
tion regime we obtain the form p| 



Wab^ 



sad. 



det$. 



-, 1/2 



IdetCd.l 



exp 



KbT 



(23) 



where Ooaad. is the down "frequency" at the saddle point 
and /i is a friction constant that determines the time 
scale (its value is fixed by a comparison with MD in the 
allowed temperature region). 

All the characteristics of the model (properties of 
the connected network and parameters in the transition 
rates) are inferred from the computed properties of the 
potential energy landscape. We use the values of the 
iV = 29 system to determine local minima properties 
(energy, curvature, stress tensor) and those of = 17 
system to determine connectivity properties (energy and 
curvature of saddle points, distances and connectivity 
among the minima). The values are extracted from the 
distribution found in the previous section in the following 
way: 
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1. we extract M energy values of the minima from 
the distribution of -/V = 29 system (we exclude the 
crystalline- like configurations); 

2. we assign to each minimum a value of curvature 
Ca = det extracted from a bivariate distribution, 
thanks to the cross-correlation between energy and 
curvature; a stress tensor value is also extracted for 
each minimum; 

3. for each minimum we randomly (in the analysis of 
the energy landscape we have found no correlation 
between energies and distances among minima) ex- 
tract 20 minima connected to it, as obtained on an 
average for the system TV = 17; 

4. we define a connection matrix Kab^ containing the 
minimum steps, i.e. the number of minima crossed, 
necessary to go from a to fe; the distance matrix dab 
is Kab times the value extracted from the distribu- 
tion of the distances among connected minima for 
N= 17; 

5. for each pair of directly connected minima we de- 
termine the energy barriers A<i>f,ar. from the value 
of the distance dab'- ^^bar. — A d" {A ~ 10^ and 
a ~ 3.7, as determined for = 17 system. Fig. 5); 

6. we assign a curvature value Csad. — \ det<i>g^^ | and 
a down "frequency" uJsad. to each saddle point, from 
bivariate distribution. 

We obtain in this way a set of parameters that describe 
the model. In order to have a good statistic description 
we considered different extractions of the parameters and 
the measured quantities were obtained by averaging over 
the extractions. 



A. Test 

Before studying the dynamical properties of the model, 
we concentrate on the static behavior obtained as asymp- 
totic solution of the master equation. In this static 
regime we can determine, in a statistical mechanical ap- 
proach, the configurational partition function 



d^'^r exp[-/3$(ri,. 



,'rN] 



(24) 



By using the approximation based on the hypothesis of 
short time local harmonic vibrations around a minimum, 
and long time collective jumps among different minima, 
we obtain 



Z(/3)^ V4''-"-)(/3) exp[-/3$,] , 



(25) 



where a labels the minima and ^a^"'™''^ is the contri- 
bution of harmonic vibrations around minimum a. This 



form of the configurational partition function emerges in 
the model as the exact infinite-time limit. The harmonic 
term is easy to calculate, being a 3A^-dimensional Gaus- 
sian integral: 



dr exp 



/3 



r 



(2^)3^/2 (deta>: 



,-1/2 



(26) 



where r = (ri, ...,r37v), and r <i>jj r = Hi.m {^a)irrL fm- 
We then obtain the approximated partition function as: 



Z(/3) ^ c /3- 



J2 (det <y exp(-/3<i>,) , (27) 



from which the thermodynamical quantities can be de- 
rived, for example for the energy: E{P) — —dp In Z. To 
check the reliability of the model we compare the quan- 
tities calculated from (|2^) with those obtained through 
MD computation. In Fig. 6 we show the potential energy 
as obtained from the model (lines) and from MD (circles). 
The MD data are obtained in the following way. Start- 
ing from high temperature we quench rapidly the system 
to low temperature, entering in a glassy state; we then 
increase the temperature up to liquid phase (open cir- 
cles). The system is subsequently slowly cooled, entering 
in the supercooled regime (100 — 70 K) and obtaining at 
the end the crystal through a first order transition (full 
circles) . The lines represent the energies determined from 
the model by taking into account all the minima (dotted 
line) and only the glassy ones (full line). A good quanti- 
tative agreement is obtained between MD and the model 
as far as the temperature is lower than about 150 K, a 
temperature in the liquid phase well above the melting 
point {Tm ^ 80 K). This result supports the correctness 
of the approximation of local- vibration / collective-jumps 
in the description of a glass-former at not-too-high tem- 
perature. This static test is a good starting point to 
extend the analysis to the dynamical regime. 



B. Equilibrium properties 

We now determine the dynamical equilibrium proper- 
ties of the model. We denote with O{r{t),r{0)) a generic 
observable which depends on collective coordinates r at 
time t and at initial time t = 0. We define the statistical 
average value of O in the model as 



<0{t) >=J2Pb Y.^-bPa{t;b,0) 



(28) 



where Oab is the value of O evaluated at the minimum 
configurations a and b: Oab — C(ziajl!b)- terms of the 
eigenvalues and eigenvectors of the transition matrix W 
we have: 
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< 0{t) >= J2 exp(A„t) J2 4"^ a 



(„) 

b 



(29) 



or in terms of the eigenvectors of the symmetric matrix 
w: 

< 0{t) >= J2 exp(A„t) J2 {p'a Pt)'^' ei") el") . 



(30) 

In the following we report a detailed analysis of the equi- 
librium dynamics for a network of 400 minima, averaging 
over 50 different extractions of the parameters that define 
the model. We measure the time autocorrelation func- 
tion of the stress tensor, the shear viscosity, the structural 
relaxation times and the mass diffusion coefficient. 

We first determine the time autocorrelation functions 
of a structural quantity which is well defined in all min- 
imum configurations, i.e. the off-diagonal microscopic 
stress tensor (Q). The correlation function is 



1 



at) - -[< a^-^it) a-(0) > +(xy) + (yz)] 
The quantity Oat hi Eq. (|2|) is in this case 



Oab 



' r ^ „zx 



'{xy) + {yz)] 



(31) 



(32) 



We have measured the correlation functions for different 
temperatures, from T = 150 K to T = 20 K. In Fig. 7 we 
report the normalized correlation functions C{t)/C{Q) at 
different temperatures (open symbols) together with the 
best stretched exponential fit (lines): 



C{t) = C(0) exp [-{t/r)'' 



(33) 



Contrary to the MD computations, which result in a two- 
step behavior for the relaxation processes (one associated 
to fast local dynamics and the other to structural slow dy- 
namics, the so called a structural processes), the model 
gives only one relaxation step, associated to the struc- 
tural processes, because the model can only describes 
the long time behavior. The results we obtain with the 
present model are consistent with those of MD in the al- 
lowed region ( i.e. above T ~ 90 K, in order to avoid 
crystallization in the MD computation). In inset (a) we 
show the temperature dependence of the stretching pa- 
rameter j3k- It emerges that: the structural relaxation 
dynamics is well represented by a stretched exponential 
decay, and that the stretching parameter (3k is strongly 
temperature dependent. 

Both results are well supported by experimental |^] 
and numerical [ p^ observations. In our case the stretch- 
ing parameters [3k decreases from a value of ^ 1 at high 
temperature to 13k ~ 0.35 at low temperature, in agree- 
ment with experimental findings ||]. 

From the behavior of the correlation function C{t) we 
can obtain information about structural relaxation time 



T. The values of r (inset (6)) are determined from the 
stretched exponential fits. We obtain an increase of many 
orders of magnitude in a small temperature range, as usu- 
ally found in many glass-forming liquids. However, we do 
not find the dramatic increase of the Vogel-Tammann- 
Fulcher type expected for fragile glass- former |jll|. It is 
possible that the observed Arrhenius behavior emerges 
as a peculiar property of the model, which would mean 
that the model is unable to capture the phenomenology 
of "fragility". It is, however, possible that the Arrhenius 
law is a genuine property of glass-forming liquids with 
Lennard- Jones interaction, as supported by a compari- 
son with a MD computation in the allowed temperature 
range. It would be very interesting to compare the be- 
havior obtained from the model to the "true" (in the 
sense of MD computation) behavior in the full temper- 
ature range (this is possible only for those system that 
avoid crystallization, like suitable binary mixtures). 

From the time autocorrelation functions of the off- 
diagonal microscopic stress tensor, we can determine the 
shear viscosity as 0] 



1 



77 : 



KbTV 



dt C{t) . 



(34) 



Also in this case, like for the relaxation times, we find a 
strong increase in a small temperature range, from -q ^ 
10-2 P at T = 150 K, to 7; ~ 10" P at T = 20 K. 
Also in this case we found a good agreement between the 
model and MD for T > 90 K, giving further support to 
the model. 

The last quantity we measured in the model was the 
mass diffusion coefficient. In order to find it we determine 
the mean square displacement 



1 1 ^ 

0{t) = ^ll(t) -r(0)p = _^|.-(<) _rl(0)| 



from which we obtain the diffusion coefficient: 



D — lim 



6t 



(35) 



(36) 



To evaluate D we use the quantity Oab = d'^b/-^^ where 
dab is defined in (10). We observe again a strong increase 
with an interesting behavior not simply linear in double 
logarithmic scale |^ . 

We conclude this section by analyzing the validity of 
the Stokes-Einstein relation in the model. The Stokes- 
Einstein relation describes in a rigorous way the diffusive 
motion of a macroscopic objects in a fluids and predicts 
the following relation: 



T 

Dec - . 



(37) 



The Stokes-Einstein relation also describes fairly well the 
diffusion at atomic scale in liquids at high temperatures. 
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By lowering temperature, as observed in many experi- 
ments iQ , one usually finds a breakdown of the Eq. ( 37 ) . 
We found ||^ that at high T the model satisfies asymptot- 
ically the Stokes-Einstein relation, but upon decreasing 
T we observe a breakdown of the relation and a fit over 
the lowest temperature data of the type 



D- 



jives the value 



^ ~ 0.28 



(38) 



(39) 



This value is in fairly good agreement with experimen- 
tal results found in fragile glass-formers, like o-terphenyl 



C. Off-equilibrium properties 

Although introduced to analyze the long time dynam- 
ics, the model allows an easy computation of the short 
time, off-equilibrium dynamics. One of the main proper- 
ties of glass-formers is the very strong increase of charac- 
teristic relaxation times when temperature is lowered. If 
these times become comparable to observational times, 
the system is no more able to explore the full acces- 
sible phase space and then to reach the thermal equi- 
librium. The observed quantities are characterized by 
off-equilibrium processes. In this regime the one-time 
quantities, such as energy or time correlation functions 
with fixed initial time, can no longer describe the physics 
of the system. The usual translational time invariance, 
valid in the equilibrium regime, is no more satisfied. One 
of the most interesting consequences of that is the fact 
that the fluctuation-dissipation relation does not longer 
hold 10 . We concentrate here on this property. 

Let be H the Hamiltonian of the system and O a 
generic observable dependent on microscopic variables. 
We define the two time autocorrelation function 



C(i,t„)=<0(i) Oit^)> , 



(40) 



where we suppose t > and < ... > now means a dy- 
namical average over initial conditions. We also intro- 
duce the response function to a perturbation e, which is 
coupled to the observable O and gives rise to a pertutbed 
Hamiltonian: 



H' = H + e{t) O , 
The response is defined as 

5 < 0{t) > 



R{t, tu 



(41) 



(42) 



where again t > In the equilibrium regime the time 
translational invariance implies the validity of fluctuation 
dissipation relation 



dr 



(43) 



where t = t — Introducing the integrate response 
function x 



Eq. 



takes the form 



dt' R{t,t') 



dxjC) 
dC 



(44) 



(45) 



In the off-equilibrium regime the fluctuation-dissipation 
ratio (|4^) is no more valid. It is possible however to gen- 
eralize the ratio introducing a violation factor X[t,tw)- 
The analytical study of some generalized mean field spin 
glass models shows that the function X(i, i^) de- 
pends on time only through the correlation function C: 
X{t,tw) = X[C{t,t-uj)\. Using this property we can 
write a generalized fluctuation-dissipation ratio in the 
off-equilibrium regime 



dx{C) 
dC 



I3X{C) 



(46) 



For short times r <C t^, we have X{C) = 1 and the 
system satisfies an equilibrium-like relation, even if it 
is confined in a small phase space region. For times 
T ^ the exploration of the phase space is an off- 
equilibrium process and this implies the violation of the 
equilibrium fluctuation-dissipation ratio. In this case we 
have X{C) < 1. The very interesting relationship be- 
tween off-equilibrium and equilibrium properties of some 
generalized spin glass model, suggests that in the case of 
one step replica symmetry breaking the X{C) function 
depends only on temperature 



X{C) = m{T) 



(47) 



'=0 



and m{T) is linear in T at low temperature. It has been 
recently suggested that structural glasses present a strik- 
ing similarity with the generalized spin glass model with 
one step replica symmetry breaking [ p^ (for a recent in- 
teresting review see Coluzzi [Q). We then expect that 
also for structural glasses the violation parameter would 
show a linear temperature dependence in the violation 
region X < \. Evidence of this behavior was found in re- 
cent numerical study of binary mixtures ; we analyze 
it in our model. 

Let be 0{r_{t)) a generic observable that depends on 
collective coordinates at time t. The average value of O 
in the off-equilibrium regime in the model is 
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a.b 



(48) 



where a and b label the minima and the sum over b is now 
limited to a certain subset of minima. We chose the M' 
highest energy states. Expression ( jisl ) differs from the 
equilibrium one, as the initial states are weighted with 
a constant term (corresponding to an infinite tempera- 
ture) rather than with the Gibbs-Boltzmann equilibrium 
weight. In this way we describe an instantaneous quench 
at time t — from T = cx) to a finite temperature T (the 
T dependence is as usual in the probability pa). The sum 
restricted to the M' initial states with highest energies 
{M' < M where M is the total number of minima; in 
our case M — 400 and M' — 20) allows a better descrip- 
tion of the off-equilibrium regime. We calculate the time 
correlation functions in the model as 

< 0(t) 0{t^) ^ ^ a Ob Pb{tn,]C, 0) pa{t\ b, t^) , 



(49) 



where the sum over b is still made over the M' highest 
minima. The quantity we determine is the time autocor- 
relation function of the off-diagonal microscopic stress 
tensor cr^^: 

C{t,t^) a^-(i) a^-(t„) > - < a'-{t) > < a'^it^) > 

(50) 

The response function is determined by the perturbed 
Hamiltonian 



H' = H + e{t) CT^^ , 

where the external field e is 

for t < t„ 



(51) 



(52) 



<a'^{t)) >, - < (7^-(t)) >,=o 



(56) 



We have determined the correlation functions C{t,tw) 
and the response x(^)^u') as functions of t for different 
times tw] the temperatures we analyze are in the range 
T = 100 20 K. In determining the response functions 
we have used a value of e small enough (e = 0.1) that 
the regime is linear, as verified by trying different e val- 
ues. In Fig. 8 we report the behavior of x versus /3C 
at temperatures T = 90 K and T = 45 K, respectively. 
While at T = 90 K the relation between x and /3C is to 
a good approximation linear with slope 1 on the whole 
range (full line), at lower temperature it is evident that 
after a first linear behavior with slope 1 (full line) an ap- 
proximately linear behavior with slope < 1 takes on at 
longer times (dashed line), as theoretically and numeri- 
cally expected. Moreover the slope of the second region 
decreases by lowering temperature: in Fig. 9 we show 
the slope m of the violation region versus T . At high 
temperature the value of m is nearly 1, while below a 
temperature of about 60 70 K m decreases linearly, as 
we expect in the hypothesis of one step replica symme- 
try breaking. Fig. 9 is limited to T > 40 K, as for lower 
temperatures the m values saturate to a limiting value 
and it is no more possible to extract correct information. 
This effect is probably due to the finite size of the system, 
because the sampling of the initial off-equilibrium states 
is not exhaustive [M' = 20). In the equilibrium analysis 
the finite size effects do not show up in the temperature 
range explored, as the sampling of the initial states is 
complete (Af = M). 

In conclusion from the analysis of the off-equilibrium 
properties of the model it emerges that the deviation 
from the usual fluctuation-dissipation relation, valid in 
the equilibrium regime, is in agreement with theoreti- 
cal predictions and numerical findings in simple glass- 
formers. 



for t > t,, 



The perturbation induces a change in the energies of the 
minima: 



The response function is 

5 < ct"^(0 >, 



R{t, tw) — 



(53) 



(54) 



e=0 



The < ... >e is evaluated in the presence of the pertur- 
bation e: 



< a^^(t) >,= Yl 0-^^ Pb{t^;c,0) pl{t;b,U 



M' ^ 

a.b.c 



(55) 



where p'^ is the solution of the master equation with the 
perturbing term. For small perturbation we obtain: 



IV. CONCLUSIONS 

The very rich phenomenology of the cooling process of 
glass-forming liquids, of the glass transition and of the 
glassy systems in general, has received in the last few 
years many important theoretical, numerical and exper- 
imental contributions. The present work is concerned 
with the numerical investigation of a simple model glass, 
a Lennard-Jones system of interacting particles. The 
main aim of the work was to determine the emergent 
properties of the system at the level of potential energy 
landscape. After a detailed analysis of the topological 
properties of the potential energy surface, we introduced 
a model which reproduces the long time dynamic behav- 
ior of the system. While in the usual MD investigations 
of relaxation the computational times are proportional 
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to physical times (with computational times of the order 
10^ s one obtains physical times of the order 10~^ s for 
system of size N ~ 10^), our model allows the study at 
very long physical times in short computational times. 

We studied both equilibrium and off-equilibrium prop- 
erties. The main equilibrium results we obtained are (i) 
the stretching of the relaxation dynamics, (m) the tem- 
perature dependence of the stretching parameter, and 
(iii) the breakdown of the Stokes- Einstein relation. If 
they are genuine properties of the glassy system ana- 
lyzed, they represent intriguing and interesting results 
that open fascinating questions about the glassy and su- 
percooled liquid behavior. 

Although introduced to investigate the long time dy- 
namics, the model is also able to describe in a sim- 
ple and direct way the off equilibrium dynamics. The 
emergent violation of the fluctuation-dissipation relation 
(that holds at equilibrium) is a very interesting feature 
and supports many conjectures about the analogy be- 
tween structural glasses and some spin glass model [ pT| . 
Moreover, the appearance of a critical temperature be- 
low which the violation takes place, seems to indicate the 
existence of a transition; its relation with the equilibrium 
properties is an open and exciting question. 

In conclusion, the analyzed features of the potential en- 
ergy landscape and the emergent properties of the model 
both at and off equilibrium, seem to provide a good 
description of the glassy systems. The method is very 
powerful for the investigation of the glassy properties by 
avoiding some of the main problems usually encountered 
in numerical studies, like the very long computational 
times in the low temperature regime or the presence of 
crystal states. We hope the analysis we performed may 
constitute a promising route in the investigation of glassy 
systems. 

We acknowledge B. Coluzzi, G. Monaco, F. Sciortino 
and P. Verrocchio for useful discussions, and D. Leporini 
who kept our attention on the fractional SE relation issue. 
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TABLE I. Correlation coefEcients r for the measured quan- 
tities X and y. 
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FIGURE CAPTIONS 

Fig. 1 - Distribution of ttie potential energy (in K per particle) 
of the minima for tlie system N = 29. In ttie inset tlie 
distribution of the curvature "fmin. ■ 

Fig. 2 - An example of potential energy profile along the least 
action path (•) between two minima, compared with the 
straight path (o). 

Fig. 3 - Energy distribution of the barriers along the least action 
paths among the minima. The inset shows the distribution 

of the curvature 'jsad. of the saddle points. 

Fig. 4 - Correlation between energies and curvatures of the min- 
ima. The highest curvature values correspond to low energy 
crystal minima. 

Fig. 5 - Correlation between energies of the barriers and distances 
among minima, in a log scale. The line is the best linear fit 
in log scale, correspondent to a power law fit (the slop is 
a = 3.7). 

Fig. 6 - Potential energy versus temperature as determined from 
MD and from the model. The MD (o) are obtained heating 
the glass and (•) cooling the liquid. Dotted line refers to 
the model using all the minima, while full line using only 

the glassy minima. 

Fig. 7 - Normalized autocorrelation functions of the off-diagonal 
microscopic stress tensor versus time at different temper- 
atures, obtained from the model. In the inset (a) the T 
dependence of the stretching parameter fSn and in the inset 
(6) the relaxation time versus l/T, both obtained from the 
stretched exponential fit of the autocorrelation functions. 

Fig. 8 - The integrate response x versus f3C at temperatures T = 
90 K (o) and T = 45 K (•). The full Une is the fluctuation- 
dissipation ratio, the dotted line is the best fit of the last 
points of T = 45 K. 

Fig. 9 - The slop m in the region of the violation of the fiuctuation- 
dissipation ratio versus temperature. The straight line fits 
the data in the violation region. 
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Fig. 2 - Potential energy landscape... - Angelani et al. 
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Fig. 5 - Potential energy landscape... - Angelani et al. 
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